﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.NumericalIntegration
{
    /// <summary>
    /// This class provides a bivariate numerical integrator. For integrating in x direction the algorithm uses 
    /// the double exponential algorithm. For integrating in x direction the algorithm uses the left point principal.
    /// </summary>
    [Serializable]
    public class BivariateDECombinationIntegrator : AbstractBivariateIntegrator
    {
        /// <summary>
        /// Initializes a new instance of the <see cref="BivariateDECombinationIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The bivariate integral to integrate.</param>
        public BivariateDECombinationIntegrator(BivariateIntegral integral)
            : base(integral)
        {
        }

        /// <summary>
        /// Integrates the integral in a numerical way. After running the numerical integration
        /// it should be checked if the result has a precision error.
        /// </summary>
        /// <param name="iterations">The maximum numbers of iterations.</param>
        /// <param name="precision">The minimum precision of the numerical integration.</param>
        /// <returns>The area of the specified integral.</returns>
        public override double Integrate(int iterations, double precision)
        {
            double distance = Math.Max(this.integral.A2, this.integral.B2) -
                              Math.Min(this.integral.A2, this.integral.B2);

            return this.Integrate((int) (distance * 4), iterations, precision);
        }

        /// <summary>
        /// Integrates the integral in a numerical way. After running the numerical integration
        /// it should be checked if the result has a precision error.
        /// </summary>
        /// <param name="parts">The number of parts, in which the function should be divided at the start of the algorithm.</param>
        /// <param name="iterations">The maximum numbers of iterations.</param>
        /// <param name="precision">The minimum precision of the numerical integration.</param>
        /// <returns>The area of the specified integral.</returns>
        public override double Integrate(int parts, int iterations, double precision)
        {
            return this.Integrate(parts, iterations, 5, precision);
        }

        /// <summary>
        /// Integrates the integral in a numerical way. After running the numerical integration
        /// it should be checked if the result has a precision error.
        /// </summary>
        /// <param name="parts">The number of parts, in which the function should be divided at the start of 
        /// the algorithm (y direction).</param>
        /// <param name="iterations">The maximum numbers of iterations.</param>
        /// <param name="exponentialIterations">The maximum numbers of iterations for the double exponential 
        /// transform integration (x direction).</param>
        /// <param name="precision">The minimum precision of the numerical integration.</param>
        /// <returns>The area of the specified integral.</returns>
        public double Integrate(int parts, int iterations, int exponentialIterations, double precision)
        {
            double sign1 = 1;
            double sign2 = 1;
            this.a1 = this.integral.A1;
            this.a2 = this.integral.A2;
            this.b1 = this.integral.B1;
            this.b2 = this.integral.B2;

            if (b1 < a1)
            {
                double tempuri = b1;

                b1 = a1;
                a1 = tempuri;
                sign1 = -1;
            }

            if (b2 < a2)
            {
                double tempuri = b2;

                b2 = a2;
                a2 = tempuri;
                sign2 = -1;
            }


            int i = parts;
            double newValue = 0;
            double oldValue = this.InternalIntegrateY(i, exponentialIterations);

            iterations += parts;

            while (i < iterations)
            {
                i++;
                newValue = this.InternalIntegrateY(i, exponentialIterations);
                this.relativeError = Math.Abs(newValue - oldValue);

                if (this.relativeError < precision)
                {
                    this.neededIterations = i;
                    this.precisionError = false;

                    if (sign1 == sign2)
                    {
                        return newValue;
                    }

                    return -newValue;
                }

                oldValue = newValue;
            }

            this.neededIterations = i;
            this.precisionError = true;

            if (sign1 == sign2)
            {
                return newValue;
            }

            return -newValue;
        }

        private double InternalIntegrateY(int partsY, int exponentialIterations)
        {
            double tempuriSecondDimension = 0;
            double hSecondDimension = (this.b2 - this.a2) / (double) partsY;

            for (int i = 0; i < partsY; i++)
            {
                tempuriSecondDimension += hSecondDimension *
                                          InternalIntegrateX(this.integral.Function, exponentialIterations, a1, b1,
                                                             this.a2 + i * hSecondDimension);
            }

            return tempuriSecondDimension;
        }

        private static double InternalIntegrateX(IHardBivariateRealFunction f, int n, double a, double b, double y)
        {
            int k, l;
            double z, exz, hcos, hsin, s, dxdz, x1, x2, w;

            n = n / 2;

            double h = 5.0 / n;
            double ss = 0.5 * f.SolveAt((a + b) / 2, y);
            double ss1 = 0.0;
            double ss2 = 0.0;
            const double pi = Math.PI / 4;

            for (k = -n; k < 0; k++)
            {
                z = h * (double) k;
                exz = Math.Exp(z);
                hcos = exz + 1.0 / exz;
                hsin = exz - 1.0 / exz;
                s = Math.Exp(pi * hsin);
                w = s + 1.0 / s;
                dxdz = hcos / (w * w);
                x1 = (b * s + a / s) / w;
                x2 = (a * s + b / s) / w;

                if (x1 != a && x1 != b)
                {
                    ss1 += dxdz * f.SolveAt(x1, y);
                }

                if (x2 != a && x2 != b)
                {
                    ss2 += dxdz * f.SolveAt(x2, y);
                }
            }

            for (l = 0; l < 3; l++)
            {
                for (k = -n; k < 0; k++)
                {
                    z = h * ((double) k + 0.5);
                    exz = Math.Exp(z);
                    hcos = exz + 1.0 / exz;
                    hsin = exz - 1.0 / exz;
                    s = Math.Exp(pi * hsin);
                    w = s + 1.0 / s;
                    dxdz = hcos / (w * w);
                    x1 = (b * s + a / s) / w;
                    x2 = (a * s + b / s) / w;

                    if (x1 != a && x1 != b)
                    {
                        ss1 += dxdz * f.SolveAt(x1, y);
                    }

                    if (x2 != a && x2 != b)
                    {
                        ss2 += dxdz * f.SolveAt(x2, y);
                    }
                }

                n *= 2;
                h /= 2;
            }

            return 2 * (b - a) * pi * h * (ss + ss1 + ss2);
        }
    }
}